Introduction

This Supplemetary Material reports the code to reproduce the analyses and figures carried out in the manuscript entitled “Glacier extinction homogenizes functional diversity” by Khelidj, Cerabolini, Caccianiga & Losapio, 2023, Journal of Vegetation Science.

The scope of this tutorial is to increase its reproducibility, clarity, transparency and dissemination. Data analysis was done using R version 4.1.2. This document was compiled with the `rmarkdown package, version Version 1.2.5033. This tutorial is licensed under CC BY-NC-ND 4.0, which means you are free to share, copy and redistribute this tutorial in any medium or format under the terms of attribution of appropriate credit, non-commercial purposes and no derivatives.

The citation is:

Khelidj R, Cerabolini BEL, Caccianiga M, Losapio G. 2023. Glacier extinction homogenizes functional diversity. Journal of Vegetation Science.

Software preparation

Install R (https://www.r-project.org) (R Core Team 2022) if you do not have it already. Then, install and load the following packages.

library(readxl)
library(FD)
library(ggplot2)
library(ggpubr)
library(lme4)
library(lmerTest)
library(car)
library(effects)
library(effectsize)
library(parameters)
library(merDeriv)
library(mgcv)
library(adiv)
library(maditr)

Data preparation

Download and import the following data tables

## trait data
traitsdb = read.csv("TraitsDB.csv")
rownames(traitsdb) = traitsdb$species
traitsdb$X = NULL

# species occurrence data
# Amola glacier 
amola_matb = read.csv('amola_matb.csv')
amola_env = read.csv('amola_env.csv')

# Cedec glacier 
cedec_matb = read.csv("cedec_matb.csv")
cedec_env = read.csv("cedec_envb.csv")

# Rutor glacier 
rutor_matb = read.csv("rutor_matb.csv")
rutor_matb = read.csv("rutor_env.csv")

# Trobio glacier 
trobio_matb = read.csv("trobio_matb.csv")
trobio_matb = read.csv("trobio_env.csv")

###
matb = list(amola_matb,cedec_matb,rutor_matb,trobio_matb)

Transform glacier retreat data

amola_env$lyears = log(amola_env$years+1)
cedec_env$lyears = log(cedec_env$years+1)
rutor_env$lyears = log(rutor_env$years+1)
trobio_env$lyears = log(trobio_env$years+1)

###
matenv = list(amola_env,cedec_env,rutor_env,trobio_env)

Calculate CV

### coefficient of variation
traitsdb$LNC.se = ifelse(traitsdb$LNC.se==0, NA, traitsdb$LNC.se)
traitsdb$LCC.se = ifelse(traitsdb$LCC.se==0, NA, traitsdb$LCC.se)
traitsdb$CN.se = ifelse(traitsdb$CN.se==0, NA, traitsdb$CN.se)
traitsdb$LA.se = ifelse(traitsdb$LA.se==0, NA, traitsdb$LA.se)
traitsdb$LDW.se = ifelse(traitsdb$LDW.se==0, NA, traitsdb$LDW.se)
traitsdb$LFW.se = ifelse(traitsdb$LFW.se==0, NA, traitsdb$LFW.se)
traitsdb$SLA.se = ifelse(traitsdb$SLA.se==0, NA, traitsdb$SLA.se)
traitsdb$LDMC.se = ifelse(traitsdb$LDMC.se==0, NA, traitsdb$LDMC.se)
traitsdb$CAN.se = ifelse(traitsdb$CAN.se==0, NA, traitsdb$CAN.se)

#
traitsdb$lnc_cv = 1/traitsdb$LNC.med * sqrt(traitsdb$LNC.n) * traitsdb$LNC.se
traitsdb$lcc_cv = 1/traitsdb$LCC.med * sqrt(traitsdb$LCC.n) * traitsdb$LCC.se
traitsdb$cn_cv = 1/traitsdb$CN.med * sqrt(traitsdb$CN.n) * traitsdb$CN.se
traitsdb$la_cv = 1/traitsdb$LA.med * sqrt(traitsdb$LA.n) * traitsdb$LA.se
traitsdb$ldw_cv = 1/traitsdb$LDW.med * sqrt(traitsdb$LDW.n) * traitsdb$LDW.se
traitsdb$ldw_cv[which(traitsdb$ldw_cv==Inf)] = NA
traitsdb$lfw_cv = 1/traitsdb$LFW.med * sqrt(traitsdb$LFW.n) * traitsdb$LFW.se
traitsdb$lfw_cv[which(traitsdb$lfw_cv==Inf)] = NA
traitsdb$sla_cv = 1/traitsdb$SLA.med * sqrt(traitsdb$SLA.n) * traitsdb$SLA.se
traitsdb$ldmc.cv = 1/traitsdb$LDMC.med * sqrt(traitsdb$LDMC.n) * traitsdb$LDMC.se
traitsdb$can_cv = 1/traitsdb$CAN.med * sqrt(traitsdb$CAN.n) * traitsdb$CAN.se

Synthesis indices

We calculate and syntheses functional diversity across plant traits, considering both trait average and trait variation.

Functional evenness and functional dispersion

We do so my means of distance-based functional diversity indices that address functional evenness (FEve) and functional dispersion (FDis). We use the R function dbFD in FD package.

for(i in 1:4){
  print(i)
  ### mean ###
  
  traits_mn = traitsdb[which(rownames(traitsdb)%in%colnames(matb[[i]])),
                         c(2,3,4,5,11,14,17,20,23,26,29,32,38)]
  
  ##### observed trends in average traits
  
  fcomp = dbFD(traits_mn, matb[[i]], corr = 'none')
  
  if(i==1){
    alltr_mni = data.frame(cbind(matenv[[i]]$lyears, fcomp$FEve, fcomp$FDis, 
                                 rep(i,nrow(matenv[[i]]))))
  }
  if(i>1){
    m0 = data.frame(cbind(matenv[[i]]$lyears, fcomp$FEve, fcomp$FDis, 
                          rep(i,nrow(matenv[[i]]))))
    alltr_mni = data.frame(rbind(alltr_mni,m0))
  }
  
  ### cv ###
  
  traits_cv = traitsdb[which(rownames(traitsdb)%in%colnames(matb[[i]])),
                         41:49]
  
  ##### observed trends in trait cv
  
  fcomp = dbFD(traits_cv, matb[[i]], corr = 'none')
  
  if(i==1){
    alltr_cvi = data.frame(cbind(matenv[[i]]$lyears, fcomp$FEve, fcomp$FDis, 
                                 rep(i,nrow(matenv[[i]]))))
  }
  if(i>1){
    m0 = data.frame(cbind(matenv[[i]]$lyears, fcomp$FEve, fcomp$FDis, 
                          rep(i,nrow(matenv[[i]]))))
    alltr_cvi = data.frame(rbind(alltr_cvi,m0))
  }
}

#
colnames(alltr_mni)[1] = 'years'
colnames(alltr_mni)[2] = 'feve'
colnames(alltr_mni)[3] = 'fdis'
colnames(alltr_mni)[4] = 'site'

alltr_mni$site = as.factor(alltr_mni$site)

alltr_mni$sprich = 0

colnames(alltr_cvi)[1] = 'years'
colnames(alltr_cvi)[2] = 'feve'
colnames(alltr_cvi)[3] = 'fdis'
colnames(alltr_cvi)[4] = 'site'

alltr_cvi$site = as.factor(alltr_cvi$site)

Evaluate relationship with species richness

alltr_cvi$sprich = 0

### relationship with alpha diversity ####
for(i in 1:4){
  print(i)
  alltr_mni$sprich[which(alltr_mni$site==i)] = rowSums(matb[[i]])
  alltr_cvi$sprich[which(alltr_cvi$site==i)] = rowSums(matb[[i]])
}

Run statistical analyses

We then look how functional evenness and functional dispersion change with glacier retreat.

Functional evenness of trait average

#### stat model ###
# mean feve
mod_feve_mni = lmer (feve ~ poly(years,2)*sprich + (1|site),
                      data = alltr_mni)

anova(mod_feve_mni)
cohens_f(mod_feve_mni, partial = TRUE)

modpars = model_parameters(mod_feve_mni, effects = "fixed", ci_method = "satterthwaite")

effsz = t_to_r(c(modpars$t[1],modpars$t[2],modpars$t[3],
                 modpars$t[4],modpars$t[5],modpars$t[6]),
       df_error = c(modpars$df[1],modpars$df[2],modpars$df[3],
                    modpars$df[4],modpars$df[5],modpars$df[6]))

Print model output and effect size

options(width = 100)
print(modpars)
## # Fixed Effects
## 
## Parameter                   | Coefficient |       SE |         95% CI |     t |     df |      p
## -----------------------------------------------------------------------------------------------
## (Intercept)                 |        0.29 |     0.03 | [ 0.21,  0.37] |  8.49 |   7.74 | < .001
## years [1st degree]          |       -1.47 |     0.24 | [-1.95, -0.99] | -6.10 | 139.96 | < .001
## years [2nd degree]          |        1.47 |     0.18 | [ 1.12,  1.82] |  8.35 | 141.39 | < .001
## sprich                      |   -9.79e-03 | 1.51e-03 | [-0.01, -0.01] | -6.48 | 125.24 | < .001
## years [1st degree] * sprich |        0.10 |     0.02 | [ 0.06,  0.13] |  5.68 | 140.14 | < .001
## years [2nd degree] * sprich |       -0.09 |     0.02 | [-0.12, -0.06] | -5.93 | 140.89 | < .001
## 
## Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald
##   t-distribution with Satterthwaite approximation.
print(effsz)
## r     |         95% CI
## ----------------------
## 0.95  | [ 0.81,  0.98]
## -0.46 | [-0.57, -0.32]
## 0.57  | [ 0.46,  0.66]
## -0.50 | [-0.61, -0.36]
## 0.43  | [ 0.29,  0.55]
## -0.45 | [-0.56, -0.31]

Functional evenness of trait variation

# cv feve
mod_feve_cvi = lmer (feve ~ poly(years,2)*sprich + (1|site),
                     data = alltr_cvi)

anova(mod_feve_cvi)
cohens_f(mod_feve_cvi, partial = TRUE)

modpars = model_parameters(mod_feve_cvi,
                           effects = "fixed",
                           ci_method = "satterthwaite")

effsz = t_to_r(c(modpars$t[1],modpars$t[2],modpars$t[3],
                 modpars$t[4],modpars$t[5],modpars$t[6]),
               df_error = c(modpars$df[1],modpars$df[2],modpars$df[3],
                            modpars$df[4],modpars$df[5],modpars$df[6]))

Print model output and effect size

print(modpars)
## # Fixed Effects
## 
## Parameter                   | Coefficient |       SE |         95% CI |     t |     df |      p
## -----------------------------------------------------------------------------------------------
## (Intercept)                 |        0.29 |     0.03 | [ 0.21,  0.37] |  8.49 |   7.74 | < .001
## years [1st degree]          |       -1.47 |     0.24 | [-1.95, -0.99] | -6.10 | 139.96 | < .001
## years [2nd degree]          |        1.47 |     0.18 | [ 1.12,  1.82] |  8.35 | 141.39 | < .001
## sprich                      |   -9.79e-03 | 1.51e-03 | [-0.01, -0.01] | -6.48 | 125.24 | < .001
## years [1st degree] * sprich |        0.10 |     0.02 | [ 0.06,  0.13] |  5.68 | 140.14 | < .001
## years [2nd degree] * sprich |       -0.09 |     0.02 | [-0.12, -0.06] | -5.93 | 140.89 | < .001
## 
## Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald
##   t-distribution with Satterthwaite approximation.
print(effsz)
## r     |         95% CI
## ----------------------
## 0.95  | [ 0.81,  0.98]
## -0.46 | [-0.57, -0.32]
## 0.57  | [ 0.46,  0.66]
## -0.50 | [-0.61, -0.36]
## 0.43  | [ 0.29,  0.55]
## -0.45 | [-0.56, -0.31]

Functional dispersion of trait average

# mean fdis
mod_fdis_mni = lmer (fdis ~ poly(years,2)*sprich + (1|site),
                     data = alltr_mni)

anova(mod_fdis_mni)
cohens_f(mod_fdis_mni, partial = TRUE)

modpars = model_parameters(mod_fdis_mni,
                           effects = "fixed",
                           ci_method = "satterthwaite")

effsz = t_to_r(c(modpars$t[1],modpars$t[2],modpars$t[3],
                 modpars$t[4],modpars$t[5],modpars$t[6]),
               df_error = c(modpars$df[1],modpars$df[2],modpars$df[3],
                            modpars$df[4],modpars$df[5],modpars$df[6]))

Print model output and effect size

options(width = 100)
print(modpars)
## # Fixed Effects
## 
## Parameter                   | Coefficient |       SE |         95% CI |     t |     df |      p
## -----------------------------------------------------------------------------------------------
## (Intercept)                 |        0.29 |     0.03 | [ 0.21,  0.37] |  8.49 |   7.74 | < .001
## years [1st degree]          |       -1.47 |     0.24 | [-1.95, -0.99] | -6.10 | 139.96 | < .001
## years [2nd degree]          |        1.47 |     0.18 | [ 1.12,  1.82] |  8.35 | 141.39 | < .001
## sprich                      |   -9.79e-03 | 1.51e-03 | [-0.01, -0.01] | -6.48 | 125.24 | < .001
## years [1st degree] * sprich |        0.10 |     0.02 | [ 0.06,  0.13] |  5.68 | 140.14 | < .001
## years [2nd degree] * sprich |       -0.09 |     0.02 | [-0.12, -0.06] | -5.93 | 140.89 | < .001
## 
## Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald
##   t-distribution with Satterthwaite approximation.
print(effsz)
## r     |         95% CI
## ----------------------
## 0.95  | [ 0.81,  0.98]
## -0.46 | [-0.57, -0.32]
## 0.57  | [ 0.46,  0.66]
## -0.50 | [-0.61, -0.36]
## 0.43  | [ 0.29,  0.55]
## -0.45 | [-0.56, -0.31]

Functional dispersion of trait variation

# cv fdis
mod_fdis_cvi = lmer (fdis ~ poly(years,2)*sprich + (1|site),
                     data = alltr_cvi)

anova(mod_fdis_cvi)
cohens_f(mod_fdis_cvi, partial = TRUE)

modpars = model_parameters(mod_fdis_cvi,
                           effects = "fixed",
                           ci_method = "satterthwaite")

effsz = t_to_r(c(modpars$t[1],modpars$t[2],modpars$t[3],
                 modpars$t[4],modpars$t[5],modpars$t[6]),
               df_error = c(modpars$df[1],modpars$df[2],modpars$df[3],
                            modpars$df[4],modpars$df[5],modpars$df[6]))

Print model output and effect size

print(modpars)
## # Fixed Effects
## 
## Parameter                   | Coefficient |       SE |         95% CI |     t |     df |      p
## -----------------------------------------------------------------------------------------------
## (Intercept)                 |        0.29 |     0.03 | [ 0.21,  0.37] |  8.49 |   7.74 | < .001
## years [1st degree]          |       -1.47 |     0.24 | [-1.95, -0.99] | -6.10 | 139.96 | < .001
## years [2nd degree]          |        1.47 |     0.18 | [ 1.12,  1.82] |  8.35 | 141.39 | < .001
## sprich                      |   -9.79e-03 | 1.51e-03 | [-0.01, -0.01] | -6.48 | 125.24 | < .001
## years [1st degree] * sprich |        0.10 |     0.02 | [ 0.06,  0.13] |  5.68 | 140.14 | < .001
## years [2nd degree] * sprich |       -0.09 |     0.02 | [-0.12, -0.06] | -5.93 | 140.89 | < .001
## 
## Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald
##   t-distribution with Satterthwaite approximation.
print(effsz)
## r     |         95% CI
## ----------------------
## 0.95  | [ 0.81,  0.98]
## -0.46 | [-0.57, -0.32]
## 0.57  | [ 0.46,  0.66]
## -0.50 | [-0.61, -0.36]
## 0.43  | [ 0.29,  0.55]
## -0.45 | [-0.56, -0.31]

Make figures

##### plots #####
plot(Effect(c("years"), mod_feve_mni))

plot(Effect(c("years"), mod_fdis_mni))

plot(Effect(c("years"), mod_feve_cvi))

plot(Effect(c("years"), mod_fdis_cvi))

plot(Effect(c("sprich"), mod_feve_mni))

plot(Effect(c("sprich"), mod_fdis_mni))

plot(Effect(c("sprich"), mod_feve_cvi))

plot(Effect(c("sprich"), mod_fdis_cvi))

plot(Effect(c("sprich", "years"), mod_feve_mni))

plot(Effect(c("sprich", "years"), mod_fdis_mni))

plot(Effect(c("sprich", "years"), mod_feve_cvi))

plot(Effect(c("sprich", "years"), mod_fdis_cvi))

### gg plots
a1 = ggplot(alltr_mni, aes(x = years,
                           y = feve,
                           colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
a2 = ggplot(alltr_mni, aes(x = years,
                           y = fdis,
                           colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
a3 = ggplot(alltr_cvi, aes(x = years,
                           y = feve,
                           colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
a4 = ggplot(alltr_cvi, aes(x = years,
                           y = fdis,
                           colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()

ggarrange(a1, a3, a2, a4,
          ncol = 2, nrow = 2)
## Warning: Removed 10 rows containing non-finite values (stat_smooth).
## Warning: Removed 10 rows containing missing values (geom_point).
## Warning: Removed 10 rows containing non-finite values (stat_smooth).
## Warning: Removed 10 rows containing missing values (geom_point).

### sprich - fdiv

##### plots #####
a1 = ggplot(alltr_mni, aes(x = sprich,
                           y = feve,
                           colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
a2 = ggplot(alltr_mni, aes(x = sprich,
                           y = fdis,
                           colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
a3 = ggplot(alltr_cvi, aes(x = sprich,
                           y = feve,
                           colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
a4 = ggplot(alltr_cvi, aes(x = sprich,
                           y = fdis,
                           colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()

ggarrange(a1, a3, a2, a4,
          ncol = 2, nrow = 2)
## Warning: Removed 10 rows containing non-finite values (stat_smooth).
## Removed 10 rows containing missing values (geom_point).
## Warning: Removed 10 rows containing non-finite values (stat_smooth).
## Warning: Removed 10 rows containing missing values (geom_point).

Functional redundancy and uniquenness

We then consider functional redundancy and uniquenness for each glacier foreland site.

#### functional redundancy ########

## amola
amola_distm = gowdis(traits_mn[which(rownames(traits_mn)%in%colnames(amola_matb)),])

amola_uniq = uniqueness(amola_matb[,which(colnames(amola_matb)%in%rownames(traits_mn))],
                        amola_distm,
                        abundance = F)

str(amola_uniq)

amola_env = cbind(amola_env,amola_uniq$red[,c(4,6)])

## cedec
cedec_distm = gowdis(traits_mn[which(rownames(traits_mn)%in%colnames(cedec_matb)),])

cedec_uniq = uniqueness(cedec_matb[,which(colnames(cedec_matb)%in%rownames(traits_mn))],
                        cedec_distm,
                        abundance = F)

cedec_uniq$Kbar

cedec_env = cbind(cedec_env,cedec_uniq$red[,c(4,6)])

## rutor
rutor_distm = gowdis(traits_mn[which(rownames(traits_mn)%in%colnames(rutor_matb)),])

rutor_uniq = uniqueness(rutor_matb[,which(colnames(rutor_matb)%in%rownames(traits_mn))],
                        rutor_distm,
                        abundance = F)

rutor_uniq$Kbar

rutor_env = cbind(rutor_env,rutor_uniq$red[,c(4,6)])

## trobio
trobio_distm = gowdis(traits_mn[which(rownames(traits_mn)%in%colnames(trobio_matb)),])

trobio_uniq = uniqueness(trobio_matb[,which(colnames(trobio_matb)%in%rownames(traits_mn))],
                         trobio_distm,
                        abundance = F)

trobio_uniq$Kbar

trobio_env = cbind(trobio_env,trobio_uniq$red[,c(4,6)])

### pool data

alltr_mni$ustar = c(amola_env$Ustar, cedec_env$Ustar, rutor_env$Ustar, trobio_env$Ustar)

Statistical analysis

mod_ustar = lmer (ustar ~ poly(years,2)*sprich + (1|site),
                     data = alltr_mni)

anova(mod_ustar)
cohens_f(mod_ustar, partial = TRUE)

modpars = model_parameters(mod_ustar, effects = "fixed", ci_method = "satterthwaite")

effsz = t_to_r(c(modpars$t[1],modpars$t[2],modpars$t[3],
                 modpars$t[4],modpars$t[5],modpars$t[6]),
               df_error = c(modpars$df[1],modpars$df[2],modpars$df[3],
                            modpars$df[4],modpars$df[5],modpars$df[6]))

Print model output

print(modpars)
## # Fixed Effects
## 
## Parameter                   | Coefficient |       SE |         95% CI |     t |     df |      p
## -----------------------------------------------------------------------------------------------
## (Intercept)                 |        0.29 |     0.03 | [ 0.21,  0.37] |  8.49 |   7.74 | < .001
## years [1st degree]          |       -1.47 |     0.24 | [-1.95, -0.99] | -6.10 | 139.96 | < .001
## years [2nd degree]          |        1.47 |     0.18 | [ 1.12,  1.82] |  8.35 | 141.39 | < .001
## sprich                      |   -9.79e-03 | 1.51e-03 | [-0.01, -0.01] | -6.48 | 125.24 | < .001
## years [1st degree] * sprich |        0.10 |     0.02 | [ 0.06,  0.13] |  5.68 | 140.14 | < .001
## years [2nd degree] * sprich |       -0.09 |     0.02 | [-0.12, -0.06] | -5.93 | 140.89 | < .001
## 
## Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald
##   t-distribution with Satterthwaite approximation.
print(effsz)
## r     |         95% CI
## ----------------------
## 0.95  | [ 0.81,  0.98]
## -0.46 | [-0.57, -0.32]
## 0.57  | [ 0.46,  0.66]
## -0.50 | [-0.61, -0.36]
## 0.43  | [ 0.29,  0.55]
## -0.45 | [-0.56, -0.31]

Figures

### Make figures
ggplot(amola_env, aes(x = lyears,
                      y = Ustar)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
## Warning: Removed 1 rows containing non-finite values (stat_smooth).
## Warning: Removed 1 rows containing missing values (geom_point).

ggplot(cedec_env, aes(x = lyears,
                      y = Ustar)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()

ggplot(rutor_env, aes(x = lyears,
                      y = Ustar)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()

ggplot(trobio_env, aes(x = lyears,
                      y = Ustar)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()

ggplot(alltr_mni, aes(x = years,
                           y = ustar,
                           colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
## Warning: Removed 1 rows containing non-finite values (stat_smooth).
## Removed 1 rows containing missing values (geom_point).

Functional diversity at the trait level

We then look at trait average and trait variation for each single trait separately.

for(i in 1:4){
  print(i)
### mean ###

traits_mn = traitsdb[which(rownames(traitsdb)%in%colnames(matb[[i]])),
                             c(3,4,5,11,14,17,20,23,26,29,32,38)]

##### observed trends in average traits

fcomp = functcomp(traits_mn, matb[[i]])

if(i==1){
alltr_mn = cbind(matenv[[i]]$lyears, fcomp, rep(i,nrow(matenv[[i]])))
}
if(i>1){
m0 = cbind(matenv[[i]]$lyears, fcomp, rep(i,nrow(matenv[[i]])))
alltr_mn = rbind(alltr_mn,m0)
}

### cv ###

traits_cv = traitsdb[which(rownames(traitsdb)%in%colnames(matb[[i]])),
                       41:49]

##### observed trends in trait cv

fcomp = functcomp(traits_cv, matb[[i]])

if(i==1){
  alltr_cv = cbind(matenv[[i]]$lyears, fcomp, rep(i,nrow(matenv[[i]])))
}
if(i>1){
  m0 = cbind(matenv[[i]]$lyears, fcomp, rep(i,nrow(matenv[[i]])))
  alltr_cv = rbind(alltr_cv,m0)
}
}

#
colnames(alltr_mn)[1] = 'years'
colnames(alltr_mn)[2] = 'latspr'
colnames(alltr_mn)[3] = 'flost'
colnames(alltr_mn)[4] = 'floper'
colnames(alltr_mn)[5] = 'lnc'
colnames(alltr_mn)[6] = 'lcc'
colnames(alltr_mn)[7] = 'cn'
colnames(alltr_mn)[8] = 'la'
colnames(alltr_mn)[9] = 'ldw'
colnames(alltr_mn)[10] = 'lfw'
colnames(alltr_mn)[11] = 'sla'
colnames(alltr_mn)[12] = 'ldmc'
colnames(alltr_mn)[13] = 'can'
colnames(alltr_mn)[14] = 'site'

alltr_mn$site = as.factor(alltr_mn$site)

#
colnames(alltr_cv)[1] = 'years'
colnames(alltr_cv)[2] = 'lnc'
colnames(alltr_cv)[3] = 'lcc'
colnames(alltr_cv)[4] = 'cn'
colnames(alltr_cv)[5] = 'la'
colnames(alltr_cv)[6] = 'ldw'
colnames(alltr_cv)[7] = 'lfw'
colnames(alltr_cv)[8] = 'sla'
colnames(alltr_cv)[9] = 'ldmc'
colnames(alltr_cv)[10] = 'can'
colnames(alltr_cv)[11] = 'site'

alltr_cv$site = as.factor(alltr_cv$site)

Statistical analysis

## mixed models & effect size####
## mean
mod.alltr.mn = rep(list(NA),12)
st.alltr.mn = data.frame(trait = gl(12,2,24, labels=colnames(alltr_mn)[2:13]),
                         estimate = NA,
                         stderr = NA,
                         bcil = NA,
                         bciu = NA,
                         df = NA,
                         tval = NA,
                         pval = NA,
                         choenr = NA,
                         rcil = NA,
                         rciu = NA)

for(i in 1:12){
  print(i)
flamod = as.formula(paste(paste(colnames(alltr_mn)[i+1], "~"),
                 paste('poly(years,2)', '(1|site)', sep="+")))

mod.alltr.mn[[i]] = lmer(flamod, data = alltr_mn)

modpars = model_parameters(mod.alltr.mn[[i]],
                           effects = "fixed",
                           ci_method = "satterthwaite")
st.alltr.mn[((i-1)*2+1),2:8] = as.numeric(modpars[2,c(2:9)])[-2]
st.alltr.mn[((i-1)*2+2),2:8] = as.numeric(modpars[3,c(2:9)])[-2]
  
effsz = t_to_r(c(modpars$t[2],modpars$t[3]),
               df_error = c(modpars$df[2],modpars$df[3]))

st.alltr.mn[((i-1)*2+1),9:11] = as.numeric(effsz[1,])[-2]
st.alltr.mn[((i-1)*2+2),9:11] = as.numeric(effsz[2,])[-2]
}

## cv
mod.alltr.cv = rep(list(NA),9)
st.alltr.cv = data.frame(trait = gl(9,2,18, labels=colnames(alltr_cv)[2:10]),
                         estimate = NA,
                         stderr = NA,
                         bcil = NA,
                         bciu = NA,
                         df = NA,
                         tval = NA,
                         pval = NA,
                         choenr = NA,
                         rcil = NA,
                         rciu = NA)

for(i in 1:9){
  print(i)
  if(i!=1  & i!=6){
  flamod = as.formula(paste(paste(colnames(alltr_cv)[i+1], "~"),
                            paste('poly(years,2)', '(1|site)', sep="+")))
  
  mod.alltr.cv[[i]] = lmer(flamod, data = alltr_cv)
  modpars = model_parameters(mod.alltr.cv[[i]],
                             effects = "fixed",
                             ci_method = "satterthwaite")
  st.alltr.cv[((i-1)*2+1),2:8] = as.numeric(modpars[2,c(2:9)])[-2]
  st.alltr.cv[((i-1)*2+2),2:8] = as.numeric(modpars[3,c(2:9)])[-2]
  
  effsz = t_to_r(c(modpars$t[2],modpars$t[3]),
                 df_error = c(modpars$df[2],modpars$df[3]))
  
  st.alltr.cv[((i-1)*2+1),9:11] = as.numeric(effsz[1,])[-2]
  st.alltr.cv[((i-1)*2+2),9:11] = as.numeric(effsz[2,])[-2]}
  if(i==1  | i==6){
    flamod = as.formula(paste(paste(colnames(alltr_cv)[i+1], "~"),
                              paste('years', '(1|site)', sep="+")))
    
    mod.alltr.cv[[i]] = lmer(flamod, data = alltr_cv)
    modpars = model_parameters(mod.alltr.cv[[i]], effects = "fixed", ci_method = "satterthwaite")
    st.alltr.cv[((i-1)*2+1),2:8] = as.numeric(modpars[2,c(2:9)])[-2]

    effsz = t_to_r(modpars$t[2],modpars$df[2])
    
    st.alltr.cv[((i-1)*2+1),9:11] = as.numeric(effsz[1,])[-2]}
}

Print model outcome

## summary table
st.alltr.mn[,-1] = round(st.alltr.mn[,-1],3)
st.alltr.cv[,-1] = round(st.alltr.cv[,-1],3)

print(st.alltr.mn)
##     trait estimate stderr     bcil    bciu     df    tval  pval choenr   rcil   rciu
## 1  latspr   -1.004   0.95   -1.540  -0.467 -3.696 144.220 0.000 -0.294 -0.429 -0.139
## 2  latspr    0.560   0.95   -0.018   1.137  1.915 145.913 0.057  0.157 -0.005  0.306
## 3   flost   -0.825   0.95   -1.149  -0.500 -5.021 143.573 0.000 -0.386 -0.507 -0.241
## 4   flost    0.005   0.95   -0.347   0.356  0.026 145.153 0.979  0.002 -0.158  0.163
## 5  floper    0.537   0.95    0.289   0.784  4.279 143.431 0.000  0.336  0.185  0.465
## 6  floper   -0.158   0.95   -0.426   0.111 -1.159 144.784 0.248 -0.096 -0.251  0.067
## 7     lnc    0.156   0.95   -0.387   0.700  0.569 144.166 0.571  0.047 -0.115  0.206
## 8     lnc   -1.284   0.95   -1.869  -0.699 -4.338 145.952 0.000 -0.338 -0.465 -0.188
## 9     lcc    9.129   0.95    7.709  10.550 12.704 144.381 0.000  0.726  0.649  0.783
## 10    lcc    3.305   0.95    1.778   4.831  4.278 145.664 0.000  0.334  0.183  0.462
## 11     cn    7.331   0.95    1.987  12.675  2.712 143.848 0.008  0.221  0.060  0.364
## 12     cn   14.132   0.95    8.361  19.903  4.840 145.824 0.000  0.372  0.226  0.494
## 13     la  561.067   0.95  237.151 884.983  3.424 143.920 0.001  0.274  0.117  0.412
## 14     la   54.999   0.95 -294.804 404.802  0.311 145.821 0.756  0.026 -0.135  0.185
## 15    ldw   37.559   0.95   19.369  55.750  4.081 143.901 0.000  0.322  0.169  0.453
## 16    ldw    2.393   0.95  -17.255  22.042  0.241 145.786 0.810  0.020 -0.141  0.179
## 17    lfw   80.421   0.95  -36.714 197.557  1.357 144.054 0.177  0.112 -0.051  0.267
## 18    lfw   16.232   0.95 -110.139 142.603  0.254 145.955 0.800  0.021 -0.140  0.180
## 19    sla   -8.929   0.95  -13.382  -4.475 -3.963 143.815 0.000 -0.314 -0.446 -0.160
## 20    sla   -9.597   0.95  -14.407  -4.786 -3.943 145.777 0.000 -0.310 -0.442 -0.157
## 21   ldmc   22.906   0.95   17.579  28.233  8.499 145.065 0.000  0.577  0.463  0.663
## 22   ldmc   11.957   0.95    6.268  17.646  4.155 142.344 0.000  0.329  0.176  0.459
## 23    can  179.062   0.95  122.597 235.527  6.268 143.591 0.000  0.464  0.329  0.572
## 24    can  -62.650   0.95 -123.768  -1.533 -2.026 145.171 0.045 -0.166 -0.315 -0.004
print(st.alltr.cv)
##    trait estimate stderr   bcil   bciu     df    tval  pval choenr   rcil   rciu
## 1    lnc    0.001   0.95  0.000  0.001  2.181 146.886 0.031  0.177  0.017  0.324
## 2    lnc       NA     NA     NA     NA     NA      NA    NA     NA     NA     NA
## 3    lcc    0.001   0.95 -0.001  0.003  1.343 143.465 0.181  0.111 -0.052  0.266
## 4    lcc   -0.005   0.95 -0.007 -0.003 -5.532 144.973 0.000 -0.418 -0.533 -0.276
## 5     cn    0.032   0.95  0.015  0.049  3.642 144.805 0.000  0.290  0.134  0.424
## 6     cn   -0.048   0.95 -0.066 -0.029 -5.087 143.215 0.000 -0.391 -0.512 -0.246
## 7     la    0.180   0.95  0.106  0.255  4.771 143.837 0.000  0.370  0.222  0.493
## 8     la   -0.216   0.95 -0.297 -0.135 -5.293 145.831 0.000 -0.401 -0.519 -0.259
## 9    ldw   -0.026   0.95 -0.126  0.073 -0.523 140.193 0.602 -0.044 -0.205  0.121
## 10   ldw    0.080   0.95 -0.039  0.198  1.326 137.155 0.187  0.112 -0.055  0.270
## 11   lfw    0.001   0.95 -0.002  0.004  0.649 146.746 0.517  0.053 -0.108  0.211
## 12   lfw       NA     NA     NA     NA     NA      NA    NA     NA     NA     NA
## 13   sla   -0.291   0.95 -0.380 -0.202 -6.465 143.966 0.000 -0.474 -0.580 -0.341
## 14   sla    0.111   0.95  0.015  0.207  2.281 145.962 0.024  0.186  0.025  0.332
## 15  ldmc   -0.249   0.95 -0.321 -0.177 -6.852 143.327 0.000 -0.497 -0.599 -0.367
## 16  ldmc    0.007   0.95 -0.071  0.085  0.180 144.446 0.857  0.015 -0.147  0.175
## 17   can    0.011   0.95 -0.078  0.100  0.245 144.236 0.807  0.020 -0.141  0.181
## 18   can    0.040   0.95 -0.055  0.135  0.832 145.809 0.407  0.069 -0.093  0.225
###

interpret_r(st.alltr.cv$choenr, rules = "gignac2016")
##  [1] "small"      NA           "small"      "large"      "moderate"   "large"      "large"     
##  [8] "large"      "very small" "small"      "very small" NA           "large"      "small"     
## [15] "large"      "very small" "very small" "very small"
## (Rules: gignac2016)

Figures

Make figures

## average ##
for(i in 1:12){
print(plot(Effect(c("years"), mod.alltr.mn[[i]])))
}

dev.off()
## null device 
##           1
####
a1 = ggplot(alltr_mn, aes(x = years,
                          y = latspr,
                          colour = site)) +
    geom_jitter(width = 0.2) +
    geom_smooth(method = lm, formula = y ~ poly(x,2)) +
    theme_classic()
a2 = ggplot(alltr_mn, aes(x = years,
                          y = flost,
                          colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
a3 = ggplot(alltr_mn, aes(x = years,
                          y = floper,
                          colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
a4 = ggplot(alltr_mn, aes(x = years,
                          y = lnc,
                          colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
a5 = ggplot(alltr_mn, aes(x = years,
                          y = lcc,
                          colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
a6 = ggplot(alltr_mn, aes(x = years,
                          y = cn,
                          colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
a7 = ggplot(alltr_mn, aes(x = years,
                          y = la,
                          colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
a8 = ggplot(alltr_mn, aes(x = years,
                          y = ldw,
                          colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
a9 = ggplot(alltr_mn, aes(x = years,
                          y = lfw,
                          colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
a10 = ggplot(alltr_mn, aes(x = years,
                          y = sla,
                          colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
a11 = ggplot(alltr_mn, aes(x = years,
                          y = ldmc,
                          colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
a12 = ggplot(alltr_mn, aes(x = years,
                          y = can,
                          colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()

ggarrange(a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12,
          ncol = 3, nrow = 4)

### cv
for(i in 1:9){
  print(plot(Effect(c("years"), mod.alltr.cv[[i]])))
}

dev.off()
## null device 
##           1
a1 = ggplot(alltr_cv, aes(x = years,
                          y = lnc,
                          colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
a2 = ggplot(alltr_cv, aes(x = years,
                          y = lcc,
                          colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
a3 = ggplot(alltr_cv, aes(x = years,
                          y = cn,
                          colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
a4 = ggplot(alltr_cv, aes(x = years,
                          y = la,
                          colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
a5 = ggplot(alltr_cv, aes(x = years,
                          y = ldw,
                          colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
a6 = ggplot(alltr_cv, aes(x = years,
                          y = lfw,
                          colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
a7 = ggplot(alltr_cv, aes(x = years,
                          y = sla,
                          colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
a8 = ggplot(alltr_cv, aes(x = years,
                          y = ldmc,
                          colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()
a9 = ggplot(alltr_cv, aes(x = years,
                          y = can,
                          colour = site)) +
  geom_jitter(width = 0.2) +
  geom_smooth(method = lm, formula = y ~ poly(x,2)) +
  theme_classic()

ggarrange(a1, a2, a3, a4, a5, a6, a7, a8, a9,
          ncol = 3, nrow = 3)
## Warning: Removed 5 rows containing non-finite values (stat_smooth).
## Warning: Removed 5 rows containing missing values (geom_point).